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Abstract 

We determine the critical layer thickness for the appearance of misfit dislocations as a 
function of the misfit e between the lattice constants of the substrate and the adsorbate 
from Kinetic Monte Carlo (KMC) simulations of heteroepitaxial growth. 
To this end, an algorithm is introduced which allows the off-lattice simulation of various 
phenomena observed in heteroepitaxial growth (see e.g. ||) including critical layer 
thickness for the appearance of misfit dislocations, or self- assembled island formation. The 
only parameters of the model are deposition flux, temperature and a pairwise interaction 
potential between the particles of the system. 

Our results are compared with a theoretical treatment of the problem and show good 
agreement with a simple power law. 

1 Introduction 

Heteroepitaxial growth has been a field of intense study in recent years. This is mainly 
due to the improved performance of semiconductor and opto-electronic devices which can be 
achieved using strained layer epitaxy. Here the heteroepitaxial growth by depositing material 
onto a substrate with the same crystal structure as the adsorbate but a slightly different 
lattice constant is of special interest. 

In the early stages of this kind of heteroepitaxial growth the adsorbate is coherent with 
the substrate. In this state the crystal topology is that of a perfect crystal, i.e. each particle 
has the same coordination number and its nearest and next-nearest neighbors form the same 
geometrical figure with only slightly modified distances M, j2|. 

As the thickness of the adsorbate film increases the elastic energy of the film rises until 
it is energetically favorable to form dislocations in order to relieve the strain. In this new in- 
coherent state the crystal topology is perturbed near the substrate/adsorbate interface. The 
thickness of the adsorbate film at which this occurs is called the critical layer thickness h c - 
Theoretical models were proposed in order to determine h c as a function of material param- 
eters fl|, ||, H|, [|]. In this letter we determine the critical layer thickness as a function of 
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the misfit between the lattice constants of the substrate and the adsorbate using computer 
simulations. 

In order to simulate heteroepitaxial growth one has to overcome the limitations of a fixed 
lattice. One obvious way to perform off-lattice simulations are molecular dynamics (MD), 
see for instance J3J] in the context of misfit dislocations. MD methods have the advantage 
of conceptual simplicity but can only be applied to rather small system sizes at very high 
temperature and deposition rates. Because of the high computational efforts the critical layer 
thickness has been determined for only a few values of the misfit. 

Here, we propose a KMC algorithm for the simulation of the early stages of heteroepitaxial 
growth. In contrast to similar off-lattice algorithms suggested before - for example by Faux et 
al. g § , Plotz et al. 1 or Schindler § - we are able to simulate heteroepitaxial growth for 
rather thick adsorbate layers and over a wide range of the misfit between the lattice constants 
of the substrate and the adsorbate. 

In the following we consider only growth in 1 + 1 dimensions. However, the method can be 
extended to 2 + 1 dimensional growth. To our knowledge this is the first time that the critical 
layer thickness for the appearance of dislocations is observed in Monte Carlo simulations. We 
find that our results fit well to a power law, as it has been observed for several semiconductor 
compounds ||. 



2 Method 



The aim here is to gain general insight into relevant mechanisms of heteroepitaxial growth. 
We therefore choose a simple Lennard-Jones potential 
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U ij {a) = 4U Q 
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as particle interaction, which is numerically easy to handle and saves computer time compared 
to more realistic empirical potentials like the EAM approximation (see e.g. p0[|). However, 
we focus on the observation of effects which should not depend on the particular choice of the 
potential. 

The equilibrium distance r$ between two isolated particles interacting via U^ becomes 
ro = -\/2a and is slightly smaller in the bulk material. Because of the isotropy of the Lennard- 
Jones potential the particles arrange in a triangular lattice. In order to save computer time 
the interaction potential Uij is cut off at a distance ry > 3ro. The interaction strength at 
this distance is less than 1% of the value at the equilibrium distance and can therefore be 
neglected. The interaction of two substrate particles is given by Uij (as). Two adsorbate 
particles interact via Uij (o~a) whereas we assume that a substrate and an adsorbate particle 
interact via \ (Uij (erg) + Uj (a a))- 

We would like to stress that growth is not simulated on a fixed lattice but rather two 
particles i and j are separated by a continuous distance rij. There are two possible events 
in our simulations: deposition and diffusion of adsorbate particles. The two-dimensional 
simulation cell is open in vertical and has periodic boundary conditions in lateral direction. 
Adsorbate particles are randomly deposited on the crystal surface with a rate Rj = LF, where 
L is the system size and F = ls~ 1 is the deposition flux. The rate Ri for a diffusion event i 
is given by an Arrhenius law 



E 



a.i 



Ri = v Q e k B T, ( 2 ) 
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where uq = 10 12 s _1 , E a i, T are the attempt frequency, the activation barrier for the diffusion 
step i and the simulation temperature, respectively. 

E a j is given by E a ^ = Et^ — E^, where En denotes the energy of the particle at the transition 
state and the energy at the binding state. Both are calculated for a frozen crystal using 
Brent's method [jDJ]. Here the saddle point search for the calculation of E t ^ can be replaced 
by a simpler maximum search in 1 + 1 dimensions. 

To consider the elastic deformation of the crystal after each microscopic event (diffusion or 
deposition) the total potential energy of the n particles - system 

n n 

Em = E E UiJ ( 3 ) 

i=l j=i+l 



is minimized using a conjugate gradient method [11| under variation of the coordinates of all 
particles (substrate and adsorbate) within a circle of radius 3ro around the particle where the 
event took place. In order to avoid strain caused by this local relaxation of the crystal, after 
a distinct number of microscopic events - depending on the misfit e - a minimization of Etot 
under variation of all particle coordinates is performed. 

The obtained rates for deposition and diffusion of adsorbate particles are used in a rejection 



free KMC simulation. Using a binary tree structure [12] an event i is chosen with the correct 
probability Ri/R, where 

R = R d + J2Ri (4) 

i 

is the total rate of all microscopic processes. Then this event is performed and the rates of 
all affected events are updated. Unlike in standard Monte Carlo simulations time does not 
advance linearly in discrete time steps At. Instead the time interval r between two microscopic 
processes is randomly drawn from a Poisson distribution P(t) = Re~ Rr by r = — where 
p is a uniformly distributed random number between and 1. 

Each simulation run starts with six atomic layers of substrate with a fixed bottom layer. 
The system size L (number of particles in the substrate's upper layer) is between L = 100 
and L = 200. Within this range we found no significant dependence of the results on L. 
Measuring lengths in units of as, a a is chosen between 0.85 and 1.11, so we can simulate 
heteroepitaxial growth for misfits 

£ = (5) 

between -15% and +11%. 



3 Results and Discussion 

In the following the temperature is set to T = 0.03^, which is sufficiently low compared to 
the melting temperature T m = 0.415^ of a pure two-dimensional Lennard-Jones solid ju|]. 

Heteroepitaxial growth is now simulated in order to determine the critical layer thickness 
h c for the appearance of dislocations as a function of the misfit e. 

To this end between 5 and 10 independent simulation runs are carried out for each value of e. 
We find that in each simulation run several dislocations appear simultaneously. Then, after 
the deposition of a few monolayers of adsorbate after the first appearance of a dislocations 
in the crystal the number of dislocations remains constant. The thickness of the adsorbate 
layer at which dislocations first appear is registered as h c . 
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Figure 1: Typical sections of crystals obtained in our simulations. The six bottom layers are 
the given substrate. The dislocations are marked with arrows. Left panel: perfect dislocation 
for e = +10%. Right panel: partial dislocation for e = +6%. The grey level for a particle 
indicates the particle's average distance to its nearest neighbors of the same kind: the lighter 
its grey level the more is this particle under compression. 



To prove the existence of a dislocation we determine the coordination number n c of each 
particle by calculating the Voronoy polyhedra 14, 15]. Voronoy polyhedra are a generalization 
of the Wigner-Seitz cell to a system without a fixed lattice. The number of sides of a Voronoy 
polyhedron gives the coordination number n c (n c = 6 for a particle in a perfect triangular 
lattice). A Burgers circuit [16] is drawn around regions of the crystal with n c ^ 6. A non- 
vanishing Burgers vector then indicates the appearance of a dislocation. 

Figure [l] shows sections of two crystals obtained in our simulations for left panel e = +10% 
and right panel e = +6%. The grey level for a particle in these pictures is obtained from the 
particle's average distance to its nearest neighbors of the same kind. The lighter the grey 
level the more is this particle under compression. 



3.1 Number of Dislocations 

Figure |2| shows the number of dislocations per unit length nrj/L counted for each value of 
e about 6 monolayers deposited adsorbate after the first appearance of a dislocation in the 
crystal (the maximum number of dislocations should be reached for this thickness of the ad- 
sorbate layer). The dashed line gives the theoretical number of dislocations in a system of size 
L under the assumption that n£> = L\e\ perfect dislocations can appear. Perfect dislocations 
(fig. [I] left panel) are those for which the crystal topology far from the substrate/adsorbate 
interface is the same as in the coherent state and the Burgers vector is therefore an integer 
multiple of the lattice vector. The formation of partial dislocations (fig. [l| right panel) - 
characterized by a Burgers vector which is a rational fraction of a lattice vector - causes the 
deviations from the theoretical results for —0.07 < e < —0.03 and 0.04 < e < 0.08. This is 
due to the fact, that partial dislocations are spatially more extended than perfect dislocations. 
For this reason in the case e > more and for e < less dislocations than n£> = L|e| have to 
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Figure 2: Number of dislocations per sys- 
tem size ri£)/L as a function of the misfit 
e. The error bars represent as the stan- 
dard error of the simulation results. The 
dashed line gives the theoretical number 
of perfect dislocations in a system of size 
L. 
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Figure 3: Critical thickness h c versus mis- 
fit |g| for s < (upper curve) and e > 
(lower curve) . The error bars are obtained 
as the standard error of the simulation re- 
sults. The solid lines are calculated using 
eq. (||) where a* = 0.15 for e < and 
a* = 0.05 for e > 0. 



be built when partial dislocations appear. Why partial dislocations only appear for distinct 
values of e is still unknown. 

3.2 Critical layer thickness 

Figure || shows the critical layer thickness h c plotted versus the absolute value of the misfit 
e. For —0.03 < e < 0.02 the critical thickness is too large to be observed in our simulations. 
The simulation results show a dependence of h c on the sign of the misfit. This was found 
before by L. Dong et al. §. We believe this dependence is due to the fact, that the Lennard- 
Jones potential is not harmonic. The potential is steeper in compression (e > 0) than in 
tension (e < 0), so that for e > it becomes favorable to form a dislocation for smaller values 
of h c . 

Our simulation results agree well with a power law (solid lines in fig. ||) 

h c = a*e- 3/2 (6) 

which was proposed by Cohen-Solal et al. [|| . 

There, an energy balance model is proposed for calculating the critical layer thickness in 
heteroepitaxial growth of semiconductor compounds. To this end the classical strain energy, 
without any change of the substrate or dislocation formation, and the deformation energy 
due to a full system of interfacial misfit dislocations were compared. The method yields the 
e -3/2 p Qwer j aW) w here h c depends mainly on the misfit e. Their model was compared with 
experimental data revealing an excellent agreement for IV-IV, III-V and II- VI semiconductor 
compounds with values for a* between a* = 0.15 and a* = 0.50. A nonlinear fit of our results 
yields a* = 0.15 for e < and a* = 0.05 for g > 
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4 Conclusion 



The goal of this examination was to determine the critical layer thickness for heteroepitaxial 
growth of Lennard-Jones particles. For this purpose we propose a novel KMC algorithm. 
We demonstrate, that with this algorithm it is possible to simulate the appearance of dislo- 
cations for a wide range of misfits and rather large system sizes compared to commonly used 
methods like molecular dynamics simulations. 

We find our simulation data in good agreement with a very simple e -3 / 2 power law, where 
the critical thickness depends mainly on the lattice mismatch. This law has been identified 
in various systems before || and may thus be considered as quite general. 

The developed algorithm is applicable in the simulation of various other phenomena ob- 
served in heteroepitaxial growth. In future work we will examine the 2D-3D transition in 
island growth and step-bunching on vicinal surfaces, for instance. 

We would like to thank A. Schindler for fruitful discussions about problems concerning 
the simulation of heteroepitaxial growth. FM and MA were supported by the Deutsche 
For schungsgemeinschaf t . 
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